Code: Tables and Figures

Load packages

library(tidyverse)
library(gtsummary)
library(flextable)
library(corrplot)
library(pacman)
library(devtools)
library(qgcomp)
library(forcats)
library(cowplot)
library(broom)
library(bkmr)
library(qgcompint) 
library(kohonen)
library(sommix)
library(prettydoc)

pacman::p_load(future, future.apply, reshape2, MASS, gWQS)

Descriptive Statistics

Data Exploration

Table 1

Distribution of sociodemographic characteristics and Child Behavior Checklist (CBCL) outcomes among Black mother-child pairs in metropolitan Atlanta, Georgia.

# Select variables
table_1 <- mix %>% dplyr::select(`Maternal Age`, `Maternal Education`, `Marital Status`, Parity,`Pre-pregnancy BMI`,`Substance Use`, `Insurance Type`, `Income to Poverty Ratio`, `Infant Sex`, `Mean Child Age`, Internalizing, Externalizing) 

# Set theme
theme_gtsummary_journal(journal = "lancet")

# Table
table_1 <- table_1 %>% 
  tbl_summary(type = all_dichotomous() ~ "categorical", statistic = all_continuous() ~ "{mean} ({sd})", digits = all_continuous() ~ 2) %>%
  modify_header(label = "**Demographic**")  %>%
  bold_labels()

table_1

Demographic	N = 1591
Maternal Age	24·97 (4·46)
Maternal Education
Less than some college	83 (52%)
College degree or higher	76 (48%)
Marital Status
Never married	75 (47%)
Married	84 (53%)
Parity
No prior births	65 (41%)
1 or more prior births	94 (59%)
Pre-pregnancy BMI
Underweight or normal	64 (40%)
Overweight	29 (18%)
Obese	66 (42%)
Substance Use
No	90 (57%)
Yes	69 (43%)
Insurance Type
Low-income Medicaid	69 (43%)
Pregnancy Medicaid	58 (36%)
Private	32 (20%)
Income to Poverty Ratio
<100%	76 (48%)
100% - 199%	52 (33%)
>200%	31 (19%)
Infant Sex
Male	81 (51%)
Female	78 (49%)
Mean Child Age	2·82 (0·64)
Internalizing	7·10 (5·29)
Externalizing	11·37 (7·62)
1 Mean (SD); n (%)

Table S1

Distributions of per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in serum samples obtained during the first trimester, among Black pregnant women in metropolitan Atlanta, Georgia.

# Select variables
table_s1 <- mix %>% dplyr::select(PFHXS, PFOS, PFOA, PFNA, PBDE47, PBDE99, PBDE100)
table_s1_PFAS <- atl_chemicals %>% dplyr::select(PFHXS, PFOS, PFOA, PFNA) %>% drop_na(PFHXS, PFOS, PFOA, PFNA)
table_s1_PBDE <- atl_chemicals %>% dplyr::select(PBDE47, PBDE99, PBDE100) %>% drop_na(PBDE47, PBDE99, PBDE100)

# Build function 
dist.func <- function(x) {
  c(
    geomean = round(exp(mean(log(x), na.rm = TRUE)), 2),
    geosd = round(exp(sd(log(x), na.rm = TRUE)), 2),
    min = round(min(x, na.rm = TRUE), 2),
    max = round(max(x, na.rm = TRUE), 2),
    median = round(median(x, na.rm = TRUE), 2)
  )
}

# Apply 
table_s1 <- as.data.frame(do.call(rbind, lapply(table_s1, dist.func)))
table_s1 <- table_s1 %>%
  mutate(sample = "Analytic")
table_s1_PFAS <- as.data.frame(do.call(rbind, lapply(table_s1_PFAS, dist.func)))
table_s1_PFAS <- table_s1_PFAS %>%
  mutate(sample = "Full")
table_s1_PBDE <- as.data.frame(do.call(rbind, lapply(table_s1_PBDE, dist.func)))
table_s1_PBDE <- table_s1_PBDE %>%
  mutate(sample = "Full")

table_s1 <- rbind(table_s1_PFAS, table_s1, table_s1_PBDE)
# Table
table_s1 <- table_s1 %>% 
mutate(Chemical = rownames(table_s1)) %>%
mutate("Geometric Mean (SD)" = paste0(geomean, " (", geosd, ")")) %>%
dplyr::select(Chemical, sample, `Geometric Mean (SD)`, min, median, max) %>%
  as_grouped_data("sample") %>%
  flextable() %>%
  set_header_labels(
  geomean = "Geometric Mean",
  geosd = "Geometric SD",
  min = "Minimum",
  max = "Maximum",
  median = "Median") %>% 
  font(fontname = "Tahoma") %>%
  autofit() %>%
  theme_alafoli()

table_s1

sample	Chemical	Geometric Mean (SD)	Minimum	Median	Maximum
Full
	PFHXS	1.16 (2.06)	0.08	1.24	6.17
	PFOS	1.89 (2.5)	0.01	2.12	12.42
	PFOA	0.62 (2.36)	0.01	0.70	4.42
	PFNA	0.26 (2.35)	0.01	0.30	2.27
Analytic
	PFHXS1	1.07 (2.07)	0.14	1.27	4.50
	PFOS1	2.37 (1.97)	0.14	2.53	12.42
	PFOA1	0.8 (2.12)	0.03	0.89	4.42
	PFNA1	0.33 (2.04)	0.01	0.34	2.27
	PBDE47	92.7 (1.98)	20.65	87.48	441.86
	PBDE99	23.07 (2.21)	7.81	23.26	308.29
	PBDE100	14.6 (2.85)	2.81	17.28	111.47
Full
	PBDE471	89.84 (2.08)	15.12	83.86	1,382.08
	PBDE991	23.1 (2.35)	7.81	22.11	341.98
	PBDE1001	13.32 (2.99)	2.81	15.44	289.20

Figure S1

Pearson correlation for the association between natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in first trimester serum samples, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.

# Correlation matrix and plot
mix_cor <- mix %>% dplyr::select(PFHXS, PFOS, PFOA, PFNA, PBDE47, PBDE99, PBDE100, Internalizing, Externalizing)

cor <- cor(mix_cor, method = "pearson")

figure_s3 <- corrplot(cor, method = "color", tl.col = "grey25", tl.cex = 6, cl.cex = 6, mar=c(0,0,1,0)) %>%
corrRect(c(1,5,8, 9))

Regression

Multiple Linear Regression

Model

# Total Sample
# Assign predictors, outcomes
predictors <- c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")
outcomes <- c("Internalizing", "Externalizing")

# Loop regression
models <- list()
for (i in 1:length(outcomes)) {
  for (j in 1:length(predictors)) {
    model_name <- paste0(outcomes[i], " ~ ", predictors[j])
    formula <- as.formula(paste(outcomes[i], " ~ ", predictors[j], " + `Mean Child Age`  + `Infant Sex` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use`"))
    model <- lm(formula = formula, data = mix)
    models[[model_name]] <- tidy(model, conf.int = TRUE) %>%
      mutate(outcome = outcomes[i])
  }
}
chemicals_adj <- do.call(rbind, models) 

# Format 
reg <- chemicals_adj %>% 
  filter(term %in% c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))  %>%
  dplyr::select(outcome, term, estimate, conf.low, conf.high) %>%
  mutate(sex = "All")

### By sex

# Male
## Filter for males
mix_m <- mix %>% filter(`Infant Sex` == "Male")

models <- list()
for (i in 1:length(outcomes)) {
  for (j in 1:length(predictors)) {
    model_name <- paste0(outcomes[i], " ~ ", predictors[j])
    formula <- as.formula(paste(outcomes[i], " ~ ", predictors[j], " + `Mean Child Age` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use`"))
    model <- lm(formula = formula, data = mix_m)
    models[[model_name]] <- tidy(model, conf.int = TRUE) %>%
      mutate(outcome = outcomes[i])
  }
}
chemicals_adj <- do.call(rbind, models) 

reg_male <- chemicals_adj %>% 
  filter(term %in% c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))  %>%
  dplyr::select(outcome, term, estimate, conf.low, conf.high) %>%
  mutate(sex = "Male")

# Female
## Filter for females
mix_f <- mix %>% filter(`Infant Sex` == "Female")

models <- list()
for (i in 1:length(outcomes)) {
  for (j in 1:length(predictors)) {
    model_name <- paste0(outcomes[i], " ~ ", predictors[j])
    formula <- as.formula(paste(outcomes[i], " ~ ", predictors[j], " + `Mean Child Age` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use`"))
    model <- lm(formula = formula, data = mix_f)
    models[[model_name]] <- tidy(model, conf.int = TRUE) %>%
      mutate(outcome = outcomes[i])
  }
}
chemicals_adj <- do.call(rbind, models) 

reg_female <- chemicals_adj %>% 
  filter(term %in% c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))  %>%
  dplyr::select(outcome, term, estimate, conf.low, conf.high) %>%
  mutate(sex = "Female")

Figure 1

Associations between natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in serum samples obtained during the first trimester, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia. PFAS in ng.mL and PBDEs in pg/mL.

Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use

# Format
reg$outcome <- factor(reg$outcome, levels = c("Internalizing", "Externalizing"))
reg$term <- factor(reg$term, levels = c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))
chemical_label <- c("PFHXS", "PFOS", "PFOA", "PFNA", "BDE-47", "BDE-99", "BDE-100")
names(chemical_label) <- c("PFHXS.log", "PFOS.log","PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")

# Figure
figure_1 <- ggplot(reg, aes(x = estimate, y = term, color = outcome)) +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high), size = 4, linewidth = 4) +
  geom_vline(xintercept = 0, color = "grey25", linetype = "dashed", alpha = 0.5) +  
  facet_grid(~outcome, scales = "free") +
  xlab("Beta (95% Confidence Interval)") +
  ylab("Chemical") +
  scale_x_continuous(n.breaks=4) +
  scale_y_discrete(labels = chemical_label) +
  scale_color_manual(values = c("grey25", "grey55")) +
  coord_flip() +
  theme_classic(base_family = "AppleGothic") +
  theme(title = element_text(size = 50, color = "grey25"),
        axis.title.x = element_text(size = 50, color = "grey25"),
        axis.title.y = element_text(size = 50, color = "grey25"),
        axis.text.x = element_text(size = 50, color = "grey25", angle = 45, vjust = 0.6),
        axis.text.y = element_text(size = 50, color = "grey25"),
        strip.background = element_rect(fill="white"),
        strip.text = element_text(colour = 'grey25', face = "bold", size = 50), 
        legend.position = "none", 
        strip.background.x = element_blank(),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank())

figure_1

Table S2

Linear regression associations (overall, child sex specific) between natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in serum samples obtained during the first trimester, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.

Total sample models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use. Sex specific models are adjusted for mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.

# Combine regression results
reg <- rbind(reg, reg_female, reg_male)

# Table
table_s2 <- reg %>%
  mutate_if(is.numeric, round, digits = 2) %>%
  mutate(Beta_CI = paste0(estimate, " (", conf.low, ", ", conf.high, ")"))  %>%
 mutate(term = recode(term, 
                     "PFHXS.log" = "PFHXS",
                     "PFOA.log" = "PFOA",
                     "PFNA.log" = "PFNA",
                     "PFOS.log" = "PFOS",
                     "PBDE47.log" = "BDE-47",
                     "PBDE99.log" = "BDE-99",
                     "PBDE100.log" = "BDE-100")) %>%
  dplyr::select(term, outcome, sex, Beta_CI) %>%
  pivot_wider(names_from = outcome, values_from = (Beta_CI)) %>%
  group_by(term, sex) %>%
  summarise(
    Internalizing = Internalizing,
    Externalizing = Externalizing) %>%
  as_grouped_data(groups = "term") %>%
  flextable() %>%
  add_body_row(values = c("","" ,"Beta (95% CI)", "Beta (95% CI)")) %>%
  font(fontname = "Tahoma") %>%
  set_header_labels(term = "Chemical", sex = "Sex") %>%
  autofit() %>%
  theme_alafoli()

table_s2

Chemical	Sex	Internalizing	Externalizing
		Beta (95% CI)	Beta (95% CI)
PFHXS
	All	-0.17 (-0.33, 0)	-0.21 (-0.36, -0.05)
	Female	-0.12 (-0.35, 0.11)	-0.17 (-0.37, 0.02)
	Male	-0.17 (-0.42, 0.07)	-0.22 (-0.48, 0.03)
PFOS
	All	-0.1 (-0.26, 0.06)	-0.21 (-0.36, -0.05)
	Female	-0.01 (-0.25, 0.23)	-0.13 (-0.33, 0.06)
	Male	-0.16 (-0.39, 0.06)	-0.24 (-0.48, -0.01)
PFOA
	All	-0.12 (-0.28, 0.04)	-0.17 (-0.33, -0.01)
	Female	-0.02 (-0.36, 0.32)	-0.13 (-0.42, 0.15)
	Male	-0.18 (-0.39, 0.02)	-0.21 (-0.43, 0)
PFNA
	All	-0.15 (-0.32, 0.01)	-0.21 (-0.36, -0.05)
	Female	-0.04 (-0.29, 0.2)	-0.06 (-0.27, 0.15)
	Male	-0.23 (-0.46, -0.01)	-0.31 (-0.55, -0.07)
BDE-47
	All	0.11 (-0.06, 0.28)	0.16 (-0.01, 0.32)
	Female	0.02 (-0.2, 0.25)	0 (-0.19, 0.18)
	Male	0.22 (-0.04, 0.48)	0.31 (0.04, 0.58)
BDE-99
	All	0.09 (-0.07, 0.25)	0.15 (-0.01, 0.3)
	Female	0.1 (-0.1, 0.31)	0.13 (-0.04, 0.3)
	Male	0.09 (-0.16, 0.34)	0.14 (-0.12, 0.41)
BDE-100
	All	0.07 (-0.1, 0.23)	0.15 (-0.01, 0.31)
	Female	0.03 (-0.2, 0.26)	0.05 (-0.14, 0.25)
	Male	0.08 (-0.17, 0.34)	0.17 (-0.1, 0.44)

Table S3

Linear regression associations between natural log transformed polybrominated diphenyl ethers (PBDEs), measured in serum samples obtained during the first trimester, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia. Sensitivity analysis additionally adjusted for total lipids.

# Assign predictors, outcomes
predictors <- c("PBDE47.log", "PBDE99.log", "PBDE100.log")
outcomes <- c("Internalizing", "Externalizing")

# Loop regression
models <- list()
for (i in 1:length(outcomes)) {
  for (j in 1:length(predictors)) {
    model_name <- paste0(outcomes[i], " ~ ", predictors[j])
    formula <- as.formula(paste(outcomes[i], " ~ ", predictors[j], " + `Mean Child Age`  + `Infant Sex` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use` + total_lipids"))
    model <- lm(formula = formula, data = mix)
    models[[model_name]] <- tidy(model, conf.int = TRUE) %>%
      mutate(N = nobs(model)) %>%
      mutate(outcome = outcomes[i])
  }
}
chemicals_adj_lipid <- do.call(rbind, models) 

# Format 
reg_lipid <- chemicals_adj_lipid %>% 
  filter(term %in% c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))  %>%
  dplyr::select(outcome, N, term, estimate, conf.low, conf.high)

# Combine regression results
# Table
table_s3 <- reg_lipid %>%
  mutate_if(is.numeric, round, digits = 2) %>%
  mutate(`Beta CI` = paste0(estimate, " (", conf.low, ", ", conf.high, ")"))  %>%
 mutate(term = recode(term, 
                     "PFHXS.log" = "PFHXS",
                     "PFOA.log" = "PFOA",
                     "PFNA.log" = "PFNA",
                     "PFOS.log" = "PFOS",
                     "PBDE47.log" = "BDE-47",
                     "PBDE99.log" = "BDE-99",
                     "PBDE100.log" = "BDE-100")) %>%
  dplyr::select(term, N, outcome, `Beta CI`) %>%
  pivot_wider(names_from = outcome, values_from = c("N", "Beta CI"),names_glue = "{outcome}_{.value}") %>%
  dplyr::select(term, Internalizing_N, `Internalizing_Beta CI`, Externalizing_N, `Externalizing_Beta CI`) %>%
  flextable() %>%
  separate_header() %>%
  font(fontname = "Tahoma") %>%
  set_header_labels(term = "Chemical") %>%
  autofit() %>%
  theme_alafoli()

table_s3

term	Internalizing		Externalizing
	N	Beta CI	N	Beta CI
BDE-47	154	0.12 (-0.05, 0.29)	154	0.16 (-0.01, 0.33)
BDE-99	154	0.09 (-0.07, 0.25)	154	0.14 (-0.02, 0.3)
BDE-100	154	0.09 (-0.08, 0.25)	154	0.17 (0, 0.33)

Table S4

Linear regression associations between natural log transformed per- and polyfluoroalkyl substances (PFAS) and polybrominated diphenyl ethers (PBDEs), measured in serum samples obtained during the first trimester, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia. Sensitivity analysis removing preterm births.

# Removing preterm brith 
mix_ptb <- mix %>% filter(ptb == 0)
# Assign predictors, outcomes
predictors <- c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")
outcomes <- c("Internalizing", "Externalizing")

# Loop regression
models <- list()
for (i in 1:length(outcomes)) {
  for (j in 1:length(predictors)) {
    model_name <- paste0(outcomes[i], " ~ ", predictors[j])
    formula <- as.formula(paste(outcomes[i], " ~ ", predictors[j], " + `Mean Child Age`  + `Infant Sex` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use`"))
    model <- lm(formula = formula, data = mix_ptb)
    models[[model_name]] <- tidy(model, conf.int = TRUE) %>%
      mutate(outcome = outcomes[i])
  }
}
chemicals_adj_ptb <- do.call(rbind, models) 

# Format 
reg_ptb <- chemicals_adj_ptb %>% 
  filter(term %in% c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))  %>%
  dplyr::select(outcome, term, estimate, conf.low, conf.high)

# Combine regression results
# Table
table_s4 <- reg_ptb %>%
  mutate_if(is.numeric, round, digits = 2) %>%
  mutate(Beta_CI = paste0(estimate, " (", conf.low, ", ", conf.high, ")"))  %>%
 mutate(term = recode(term, 
                     "PFHXS.log" = "PFHXS",
                     "PFOA.log" = "PFOA",
                     "PFNA.log" = "PFNA",
                     "PFOS.log" = "PFOS",
                     "PBDE47.log" = "BDE-47",
                     "PBDE99.log" = "BDE-99",
                     "PBDE100.log" = "BDE-100")) %>%
  dplyr::select(term, outcome, Beta_CI) %>%
  pivot_wider(names_from = outcome, values_from = (Beta_CI)) %>%
  group_by(term) %>%
  summarise(
    Internalizing = Internalizing,
    Externalizing = Externalizing) %>%
  flextable() %>%
  add_body_row(values = c("","Beta (95% CI)", "Beta (95% CI)")) %>%
  font(fontname = "Tahoma") %>%
  set_header_labels(term = "Chemical") %>%
  autofit() %>%
  theme_alafoli()

table_s4

Chemical	Internalizing	Externalizing
	Beta (95% CI)	Beta (95% CI)
BDE-100	0.05 (-0.13, 0.22)	0.13 (-0.04, 0.31)
BDE-47	0.17 (-0.02, 0.35)	0.16 (-0.03, 0.35)
BDE-99	0.13 (-0.04, 0.3)	0.15 (-0.03, 0.32)
PFHXS	-0.18 (-0.35, -0.01)	-0.24 (-0.41, -0.07)
PFNA	-0.16 (-0.33, 0.01)	-0.22 (-0.4, -0.05)
PFOA	-0.11 (-0.29, 0.06)	-0.2 (-0.37, -0.03)
PFOS	-0.12 (-0.29, 0.05)	-0.24 (-0.41, -0.07)

Mixture Methods

Quantile G-Computation

Total Sample

# Define exposure, outcome, dataset
mixture <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log, Internalizing, Externalizing, `Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, `Pre-pregnancy BMI`, `Substance Use`)

exposure.names<-c("PFHXS.log", "PFOS.log", "PFOA.log","PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")

outcomes <- c("Internalizing", "Externalizing")

# Build function for  qgcomp 

qgcompdata <- function(outcome) {
  model <- qgcomp.glm.noboot(as.formula(paste(outcome, "~", "PFHXS.log + PFOS.log + PFOA.log + PFNA.log + PBDE47.log + PBDE99.log + PBDE100.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` + +`Substance Use`")),
                           expnms = exposure.names,
                           dat = mixture,
                           family = gaussian(),
                           q = 4)

 coefs <- data.frame(
    estimate = round(model$psi,2),
    lower_ci = round(model$ci.coef[-1, 1],2),
    upper_ci = round(model$ci.coef[-1, 2],2),
    outcome = outcome, 
    row.names = NULL
  )
adj_df <- data.frame()
adj_df <- rbind(adj_df, coefs)
print(adj_df)
}

# Run
results_all_int <- qgcompdata(outcome = "Internalizing")
results_all_ext <- qgcompdata(outcome = "Externalizing")

# Combine, format
qgcomp_all_df <- rbind(results_all_int, results_all_ext)
qgcomp_all_df <- qgcomp_all_df %>% mutate(group = "Total Mixture")

Chemical Specific

##### PFAS

# Define exposure, outcome dataset
mixture <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log, Internalizing, Externalizing, `Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, `Pre-pregnancy BMI`, `Insurance Type`, `Substance Use`, `Marital Status`)

outcomes <- c("Internalizing", "Externalizing")

exposure.names<-c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log")

# Function for  qgcomp 

qgcompdata <- function(outcome) {
  model <- qgcomp.glm.noboot(as.formula(paste(outcome, "~", "PFHXS.log + PFOS.log + PFOA.log + PFNA.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` + `Insurance Type` +`Substance Use` +`Marital Status`")),
                  expnms = exposure.names,
                  dat = mixture,
                 family = gaussian(),
                 q = 4)

 coefs <- data.frame(
    estimate = round(model$psi,2),
    lower_ci = round(model$ci.coef[-1, 1],2),
    upper_ci = round(model$ci.coef[-1, 2],2),
    outcome = outcome, 
    row.names = NULL
  )
adj_df <- data.frame()
PFAS_df <- rbind(adj_df, coefs)
print(PFAS_df)
}

# Run
results_pfas_int <- qgcompdata(outcome = "Internalizing")
results_pfas_ext <- qgcompdata(outcome = "Externalizing")

qgcomp_pfas_df <- rbind(results_pfas_int, results_pfas_ext)

#### PBDE

# Define exposure, outcome, dataset

mixture <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log, Internalizing, Externalizing, `Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, `Pre-pregnancy BMI`, `Insurance Type`, `Substance Use`, `Marital Status`)

exposure.names<-c("PBDE47.log", "PBDE99.log", "PBDE100.log")

outcomes <- c("Internalizing", "Externalizing")

qgcompdata <- function(outcome) {
  model <- qgcomp.glm.noboot(as.formula(paste(outcome, "~", "PBDE47.log + PBDE99.log + PBDE100.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` + `Insurance Type` + +`Substance Use` +`Marital Status`")),
                  expnms = exposure.names,
                  dat = mixture,
                 family = gaussian(),
                 q = 4)

 coefs <- data.frame(
    estimate = round(model$psi,2),
    lower_ci = round(model$ci.coef[-1, 1],2),
    upper_ci = round(model$ci.coef[-1, 2],2),
    outcome = outcome, 
    row.names = NULL
  )
adj_df <- data.frame()
adj_df <- rbind(adj_df, coefs)
print(adj_df)
}

# Run
results_pbde_int <- qgcompdata(outcome = "Internalizing")
results_pbde_ext <- qgcompdata(outcome = "Externalizing")

# Combine results
qgcomp_pbde_df <- rbind(results_pbde_int, results_pbde_ext)
qgcomp_chemicals_df <- rbind(qgcomp_pbde_df, qgcomp_pfas_df)
qgcomp_chemicals_df <-  qgcomp_chemicals_df %>% mutate(group = c("PBDE", "PBDE", "PFAS", "PFAS"))

Figure 2

Quantile g-computation cumulative effect (95% confidence intervals) and weights reflecting the associations between a mixture of natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in serum samples obtained during the first trimester, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.

Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use. multiply by negative 1

##### TOTAL SAMPLE
# Forest plot
qgcomp_all_df$outcome <- factor(qgcomp_all_df$outcome, levels = c("Internalizing", "Externalizing"))
forest_all_qgcomp <- ggplot(qgcomp_all_df, aes(estimate, outcome)) +
  geom_pointrange(aes(xmin = lower_ci, xmax = upper_ci),
                  size = 2, linewidth = 2) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey25") +
  coord_flip() +
  scale_x_continuous(limits = c(-0.6, 0.4), breaks = c(-0.6, -0.4, -0.2, 0, 0.2, 0.4)) +
  ggtitle("Total Sample") + 
  theme_classic(base_family = "AppleGothic") +
  labs(x = "Beta (95 % Confidence Interval)", y = "CBCL Outcome") +
  theme(text = element_text(size = 25, color = "grey25", family = "Helvetica"),
        legend.title = element_text(color = "grey25"),
        axis.text = element_text(color = "grey25", size = 25),
        plot.title = element_text(hjust = 0.5))

# Weights
# Internalizing
exposure.names<-c("PFHXS.log", "PFOS.log", "PFOA.log","PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")

model <- qgcomp.glm.noboot(Internalizing ~ PFHXS.log + PFOS.log + PFOA.log + PFNA.log + PBDE47.log + PBDE99.log + PBDE100.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` +`Substance Use`,
                           expnms = exposure.names,
                           dat = mixture,
                           family = gaussian(),
                           q = 4)
# Extract weights
pos <- as.data.frame(model$pos.weights) 
# Format
pos <- pos %>% rename(weight = `model$pos.weights`) %>%
  mutate(chemical = c("BDE-47", "PFOS", "BDE-99", "PFOA"), outcome = "Internalizing", Direction = "Positive")
neg <- as.data.frame(model$neg.weights)
neg <- neg %>% rename(weight = `model$neg.weights`) %>%
  mutate(chemical = c("PFHXS", "PFNA", "BDE-100"), outcome = "Internalizing", Direction = "Negative")

neg$weight <- neg$weight*-1


int_all_weights_qgcomp <- rbind(pos, neg)



# Externalizing
model <- qgcomp.glm.noboot(Externalizing ~ PFHXS.log + PFOS.log + PFOA.log + PFNA.log + PBDE47.log + PBDE99.log + PBDE100.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` +`Substance Use`,
                           expnms = exposure.names,
                           dat = mixture,
                           family = gaussian(),
                           q = 4)

# Extract weights
pos <- as.data.frame(model$pos.weights) 
pos <- pos %>% rename(weight = `model$pos.weights`) %>%
  mutate(chemical = c("BDE-47", "PFOA", "BDE-100", "BDE-99", "PFOS"), outcome = "Externalizing", Direction = "Positive")
neg <- as.data.frame(model$neg.weights)
neg <- neg %>% rename(weight = `model$neg.weights`) %>%
  mutate(chemical = c("PFHXS", "PFNA"), outcome = "Externalizing", Direction = "Negative")

neg$weight <- neg$weight*-1

ext_all_weights_qgcomp <- rbind(pos, neg)

# Plot
int_all_weights_qgcomp <- ggplot(int_all_weights_qgcomp, aes(fct_reorder(chemical, weight), weight, fill = Direction)) +
  geom_col() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
  labs(x = "Chemical", y = "Weight", fill = "Direction", title = "Internalizing") +
scale_fill_manual(values = c("salmon3", "steelblue3")) +
  scale_y_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
  theme_classic(base_family = "AppleGothic") +
  theme(text = element_text(size = 25, color = "grey25"),
        legend.title = element_text(color = "grey25"),
        axis.text.x = element_text(color = "grey25", angle = 45, hjust = 1, size = 25),
        axis.text.y = element_text(color = "grey25", size = 25),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5))

ext_all_weights_qgcomp <- ggplot(ext_all_weights_qgcomp, aes(fct_reorder(chemical, weight), weight, fill = Direction)) +
  geom_col() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
  labs(x = "Chemical", y = "Weight", fill = "Direction", title = "Externalizing") +
    scale_y_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
scale_fill_manual(values = c("salmon3", "steelblue3")) +
  theme_classic(base_family = "AppleGothic") +
  theme(text = element_text(size = 25, color = "grey25"),
        legend.title = element_text(color = "grey25"),
        axis.text.x = element_text(color = "grey25", angle = 45, hjust = 1, size = 25),
        axis.text.y = element_text(color = "grey25", size = 25),
        legend.position = "bottom",
        plot.title = element_text(hjust = 0.5),
        legend.key.size = unit(3, "cm"),
        legend.text = element_text(size = 25))

##### CHEMCIAL SPECIFIC

# Forest plot
qgcomp_pfas_df$outcome <- factor(qgcomp_pfas_df$outcome, levels = c("Internalizing", "Externalizing"))
forest_pfas_qgcomp <- ggplot(qgcomp_pfas_df, aes(estimate, outcome)) +
  geom_pointrange(aes(xmin = lower_ci, xmax = upper_ci),
                  size = 2, linewidth = 2, color = "black", position = position_dodge(width = 0.25)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey25") +
  scale_x_continuous(limits = c(-0.6, 0.4), breaks = c(-0.6, -0.4, -0.2, 0, 0.2, 0.4)) +
  ggtitle("PFAS") + 
  coord_flip() +
  theme_classic(base_family = "AppleGothic") +
  labs(x = "Beta (95 % Confidence Interval)", y = "CBCL Outcome", linetype = "Chemical Class") +
  theme(text = element_text(size = 25, color = "grey25", family = "Helvetica"),
        legend.title = element_text(color = "grey25"),
        axis.text = element_text(color = "grey25", size = 25),
        plot.title = element_text(hjust = 0.5))

forest_pbde_qgcomp <- ggplot(qgcomp_pbde_df, aes(estimate, outcome)) +
  geom_pointrange(aes(xmin = lower_ci, xmax = upper_ci),
                  size = 2, linewidth = 2, color = "black", position = position_dodge(width = 0.25)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey25") +
  scale_x_continuous(limits = c(-0.6, 0.4), breaks = c(-0.6, -0.4, -0.2, 0, 0.2, 0.4)) +
  ggtitle("PBDE") +
  coord_flip() +
  theme_classic(base_family = "AppleGothic") +
  labs(x = "Beta (95 % Confidence Interval)", y = "CBCL Outcome", linetype = "Chemical Class") +
  theme(text = element_text(size = 25, color = "grey25", family = "Helvetica"),
        legend.title = element_text(color = "grey25"),
        axis.text = element_text(color = "grey25", size = 25),
        plot.title = element_text(hjust = 0.5))

# Weights
# PFAS 
exposure.names<-c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log")

model <- qgcomp.glm.noboot(Internalizing ~ PFHXS.log + PFOS.log + PFOA.log + PFNA.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` +`Substance Use`,
                           expnms = exposure.names,
                           dat = mixture,
                           family = gaussian(),
                           q = 4)
# Extract weights
pos <- as.data.frame(model$pos.weights) 
pos <- pos %>% rename(weight = `model$pos.weights`) %>%
  mutate(chemical = c("PFOS", "PFOA"), outcome = "Internalizing", Direction = "Positive")
neg <- as.data.frame(model$neg.weights)
neg <- neg %>% rename(weight = `model$neg.weights`) %>%
  mutate(chemical = c("PFNA", "PFHXS"), outcome = "Internalizing", Direction = "Negative")

neg$weight <- neg$weight*-1

weights_int_pfas_qgcomp <- rbind(pos, neg)

# Externalizing
model <- qgcomp.glm.noboot(Externalizing ~ PFHXS.log + PFOS.log + PFOA.log + PFNA.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` +`Substance Use`,
                           expnms = exposure.names,
                           dat = mixture,
                           family = gaussian(),
                           q = 4)
# Extract weights
pos <- as.data.frame(model$pos.weights) 
pos <- pos %>% rename(weight = `model$pos.weights`) %>%
  mutate(chemical = c("PFOA", "PFOS"), outcome = "Externalizing", Direction = "Positive")
neg <- as.data.frame(model$neg.weights)
neg <- neg %>% rename(weight = `model$neg.weights`) %>%
  mutate(chemical = c("PFNA", "PFHXS"), outcome = "Externalizing", Direction = "Negative")

neg$weight <- neg$weight*-1

weights_ext_pfas_qgcomp <- rbind(pos, neg)

# PBDE
exposure.names<-c("PBDE47.log", "PBDE99.log", "PBDE100.log")

model <- qgcomp.glm.noboot(Internalizing ~ PBDE47.log + PBDE99.log + PBDE100.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` +`Substance Use`,
                           expnms = exposure.names,
                           dat = mixture,
                           family = gaussian(),
                           q = 4)
pos <- as.data.frame(model$pos.weights) 
pos <- pos %>% rename(weight = `model$pos.weights`) %>%
  mutate(chemical = c("BDE-47", "BDE-99"), outcome = "Internalizing", Direction = "Positive")
neg <- as.data.frame(model$neg.weights)
neg <- neg %>% rename(weight = `model$neg.weights`) %>%
  mutate(chemical = c("BDE-100"), outcome = "Internalizing", Direction = "Negative")

neg$weight <- neg$weight*-1

# Extract weights
weights_int_pbde_qgcomp <- rbind(pos, neg)

# Externalizing
model <- qgcomp.glm.noboot(Externalizing ~ PBDE47.log + PBDE99.log + PBDE100.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` +`Substance Use`,
                           expnms = exposure.names,
                           dat = mixture,
                           family = gaussian(),
                           q = 4)
# Extract weights
pos <- as.data.frame(model$pos.weights) 
pos <- pos %>% rename(weight = `model$pos.weights`) %>%
  mutate(chemical = c("BDE-47", "BDE-99", "BDE-100"), outcome = "Externalizing", Direction = "Positive")

weights_ext_pbde_qgcomp <- pos

# Plot
int_pfas_weights_qgcomp <- ggplot(weights_int_pfas_qgcomp, aes(fct_reorder(chemical, weight), weight, fill = Direction)) +
  geom_col() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(x = "Chemical", y = "Weight", fill = "Direction", title = "Internalizing") +
scale_fill_manual(values = c("salmon3", "steelblue3")) +
  scale_y_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
  theme_classic(base_family = "AppleGothic") +
  theme(text = element_text(size = 20, color = "grey25"),
        legend.title = element_text(color = "grey25"),
        axis.text.x = element_text(color = "grey25", angle = 45, hjust = 1, size = 25),
        axis.text.y = element_text(color = "grey25", size = 25),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5))

ext_pfas_weights_qgcomp <- ggplot(weights_ext_pfas_qgcomp, aes(fct_reorder(chemical, weight), weight, fill = Direction)) +
  geom_col() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(x = "Chemical", y = "Weight", fill = "Direction", title = "Externalizing") +
scale_fill_manual(values = c("salmon3", "steelblue3")) +
    scale_y_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
  theme_classic(base_family = "AppleGothic") +
  theme(text = element_text(size = 25, color = "grey25"),
        legend.title = element_text(color = "grey25"),
        axis.text.x = element_text(color = "grey25", angle = 45, hjust = 1, size = 25),
        axis.text.y = element_text(color = "grey25", size = 25),
        legend.position = "bottom",
        plot.title = element_text(hjust = 0.5),
        legend.key.size = unit(3, "cm"),
        legend.text = element_text(size = 25))


int_pbde_weights_qgcomp <- ggplot(weights_int_pbde_qgcomp, aes(fct_reorder(chemical, weight), weight, fill = Direction)) +
  geom_col() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(x = "Chemical", y = "Weight", fill = "Direction", title = "Internalizing") +
scale_fill_manual(values = c("salmon3", "steelblue3")) +
  theme_classic(base_family = "AppleGothic") +
  scale_y_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
  theme(text = element_text(size = 25, color = "grey25"),
        legend.title = element_text(color = "grey25"),
        axis.text.x = element_text(color = "grey25", angle = 45, hjust = 1, size = 25),
        axis.text.y = element_text(color = "grey25", size = 25),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5))

ext_pbde_weights_qgcomp <- ggplot(weights_ext_pbde_qgcomp, aes(fct_reorder(chemical, weight), weight, fill = Direction)) +
  geom_col() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(x = "Chemical", y = "Weight", fill = "Direction", title = "Externalizing") +
scale_fill_manual(values = c("steelblue3")) +
  theme_classic(base_family = "AppleGothic") +
  scale_y_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
  theme(text = element_text(size = 25, color = "grey25"),
        legend.title = element_text(color = "grey25"),
        axis.text.x = element_text(color = "grey25", angle = 45, hjust = 1, size = 25),
        axis.text.y = element_text(color = "grey25", size = 25),
        legend.position = "bottom",
        plot.title = element_text(hjust = 0.5),
        legend.key.size = unit(3, "cm"),
        legend.text = element_text(size = 25))

figure_3 <- plot_grid(forest_all_qgcomp, forest_pfas_qgcomp, forest_pbde_qgcomp, int_all_weights_qgcomp,
int_pfas_weights_qgcomp, int_pbde_weights_qgcomp,ext_all_weights_qgcomp, ext_pfas_weights_qgcomp, 
ext_pbde_weights_qgcomp, ncol = 3, nrow = 3, align = "v")

figure_3

Table S5

Association between the total mixture of per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in first trimester serum samples, a mixture of PFAS alone, and a mixture of PBDEs alone, estimated using quantile g-computation, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.

Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.

# Combine results
table_s6 <- rbind(qgcomp_all_df, qgcomp_chemicals_df)

# Table
table_s6 <- table_s6 %>%
  mutate_if(is.numeric, round, digits = 2) %>%
  mutate("Beta (CI)" = paste0(estimate, " (", lower_ci, ", ", upper_ci, ")"))  %>%
  dplyr::select(group, outcome, `Beta (CI)`) %>%
  pivot_wider(names_from = outcome, values_from = `Beta (CI)`) %>%
  flextable() %>%
  font(fontname = "Tahoma") %>%
  add_body_row(values = c("","Beta (95% CI)", "Beta (95% CI)")) %>%
  merge_v(j = "group") %>%
  set_header_labels(
  group = "") %>%
  autofit() %>%
  theme_alafoli()

table_s6

	Internalizing	Externalizing
	Beta (95% CI)	Beta (95% CI)
Total Mixture	-0.04 (-0.28, 0.2)	-0.06 (-0.3, 0.18)
PBDE	0.13 (-0.03, 0.3)	0.2 (0.04, 0.36)
PFAS	-0.2 (-0.39, -0.01)	-0.28 (-0.46, -0.1)

Bayeasian Kernal Machine Regression

Total Sample

# Define dataset
mixture <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log, Internalizing, Externalizing, `Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, ppbmi, `Substance Use`)

# Define mixture, covariates
df.exp <- mixture %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log)

df.cov <- mixture %>% dplyr::select(`Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, ppbmi, `Substance Use`)

df.cov <- sapply(df.cov, as.numeric)

# BKMR


# Internalizing
set.seed(123)
fit_int <- kmbayes(y = mixture$Internalizing,
                Z = df.exp,
                X = df.cov,
                iter = 10000, 
                verbose = FALSE, 
                varsel = TRUE,
                family='gaussian')

# Extract pips
pips_int <- ExtractPIPs(fit_int)
pips_int <- pips_int %>% mutate(outcome = "Internalizing", group = "Total Sample")

# Externalizing
set.seed(123)
fit_ext <- kmbayes(y = mixture$Externalizing,
                Z = df.exp,
                X = df.cov,
                iter = 10000, 
                verbose = FALSE, 
                varsel = TRUE,
                family='gaussian')

# Extract pips
pips_ext <- ExtractPIPs(fit_ext)
pips_ext <- pips_ext %>% mutate(outcome = "Externalizing", group = "Total Sample")

Chemical Specific

## PFAS 
# Define dataset
mixture <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log, Internalizing, Externalizing, `Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, ppbmi, `Substance Use`)

# Define mixture, covariates
df.exp <- mixture %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log)

df.cov <- mixture %>% dplyr::select(`Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, ppbmi,`Substance Use`)

df.cov <- sapply(df.cov, as.numeric)

# BKMR
# Internalizing
set.seed(123)
fit_int_pfas <- kmbayes(y = mixture$Internalizing,
                Z = df.exp,
                X = df.cov,
                iter = 10000, 
                verbose = FALSE, 
                varsel = TRUE,
                family='gaussian')


# Extract pips
pips_int_pfas <- ExtractPIPs(fit_int_pfas)
pips_int_pfas <- pips_int_pfas %>% mutate(outcome = "Internalizing", group = "PFAS")


# Externalizing
set.seed(123)
fit_ext_pfas <- kmbayes(y = mixture$Externalizing,
                Z = df.exp,
                X = df.cov,
                iter = 10000, 
                verbose = FALSE, 
                varsel = TRUE,
                family='gaussian')

# Extract pips
pips_ext_pfas <- ExtractPIPs(fit_ext_pfas)
pips_ext_pfas <- pips_ext_pfas %>% mutate(outcome = "Externalizing", group = "PFAS")


## PBDE
# Define dataset
mixture <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log, Internalizing, Externalizing, `Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, ppbmi, `Substance Use`) 

# Define mixture, covariates
df.exp <- mixture %>% dplyr::select(PBDE47.log, PBDE99.log, PBDE100.log)

df.cov <- mixture %>% dplyr::select(`Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, ppbmi, `Substance Use`)

df.cov <- sapply(df.cov, as.numeric)

# BKMR


# Internalizing 
set.seed(123)
fit_int_pbde <- kmbayes(y = mixture$Internalizing,
                Z = df.exp,
                X = df.cov,
                iter = 10000, 
                verbose = FALSE, 
                varsel = TRUE,
                family='gaussian')

# Extract pips
pips_int_pbde <- ExtractPIPs(fit_int_pbde)
pips_int_pbde <- pips_int_pbde %>% mutate(outcome = "Internalizing", group = "PBDE")

# Externalizing
set.seed(123)
fit_ext_pbde <- kmbayes(y = mixture$Externalizing,
                Z = df.exp,
                X = df.cov,
                iter = 10000, 
                verbose = FALSE, 
                varsel = TRUE,
                family='gaussian')

# Extract pips
pips_ext_pbde <- ExtractPIPs(fit_ext_pbde)
pips_ext_pbde <- pips_ext_pbde %>% mutate(outcome = "Externalizing", group = "PBDE")

Table S6

Posterior inclusion probabilities (PIPs) reflecting the importance of each exposure in Bayesian Kernal Machine Regression models estimating the association between natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in first trimester serum samples, a mixture of PFAS alone, and a mixture of PBDEs alone.

Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.

# Combine results
bkmr <- rbind(pips_int, pips_ext, pips_int_pfas, pips_ext_pfas, pips_int_pbde, pips_ext_pbde)

# Table
table_s7 <- bkmr %>%
  mutate_if(is.numeric, round, digits = 2) %>%
  dplyr::select(outcome, group, variable, PIP) %>%
  mutate(variable = recode(variable, 
                     "PFHXS.log" = "PFHXS",
                     "PFOA.log" = "PFOA",
                     "PFNA.log" = "PFNA",
                     "PFOS.log" = "PFOS",
                     "PBDE47.log" = "BDE-47",
                     "PBDE99.log" = "BDE-99",
                     "PBDE100.log" = "BDE-100")) %>%
  pivot_wider(names_from = outcome, values_from = PIP) %>%
  as_grouped_data(groups = "group") %>%
  flextable() %>%
  add_body_row(values = c("", "", "PIP", "PIP")) %>%
  set_header_labels(group = "", variable = "Chemical") %>%
  font(fontname = "Tahoma") %>%
  autofit() %>%
  theme_alafoli()

table_s7

	Chemical	Internalizing	Externalizing
		PIP	PIP
Total Sample
	PFHXS	0.47	0.44
	PFOS	0.27	0.31
	PFOA	0.24	0.16
	PFNA	0.34	0.32
	BDE-47	0.46	0.28
	BDE-99	0.37	0.27
	BDE-100	0.26	0.36
PFAS
	PFHXS	0.50	0.48
	PFOS	0.22	0.32
	PFOA	0.19	0.23
	PFNA	0.35	0.41
PBDE
	BDE-47	0.54	0.44
	BDE-99	0.56	0.46
	BDE-100	0.30	0.38

Figure 3

Cumulative effect (95% credible intervals) of the BKMR mixture of natural log transformed per- and polyfluoroalkyl substances (PFAS) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in serum samples obtained during the first trimester, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.

Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.

### TOTAL SAMPLE
df.exp <- mixture %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log)

# Internalizing 
risk.overall_int <- OverallRiskSummaries(fit =fit_int, y = mixture$Internalizing, Z = as.matrix(df.exp), X = as.matrix(df.cov), 
                                      qs = seq(0.25, 0.75, by = 0.05), 
                                      q.fixed = 0.5, method = "exact")

# Externalizing

risk.overall_ext <- OverallRiskSummaries(fit = fit_ext, y = mixture$Externalizing, Z = as.matrix(df.exp), X = as.matrix(df.cov), 
                                      qs = seq(0.25, 0.75, by = 0.05), 
                                      q.fixed = 0.5, method = "exact")

# Merge
risk.overall_int <- risk.overall_int %>% mutate(outcome = "Internalizing") 
risk.overall_ext <- risk.overall_ext  %>% mutate(outcome = "Externalizing")
risk.overall <- rbind(risk.overall_int, risk.overall_ext)


## PFAS
df.exp <- mixture %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log)

# Internalizing
risk.overall_int_pfas <- OverallRiskSummaries(fit = fit_int_pfas, y = mixture$Internalizing, Z = as.matrix(df.exp), X = as.matrix(df.cov), 
                                      qs = seq(0.25, 0.75, by = 0.05), 
                                      q.fixed = 0.5, method = "exact")

# Externalizing

risk.overall_ext_pfas <- OverallRiskSummaries(fit = fit_ext_pfas, y = mixture$Externalizing, Z = as.matrix(df.exp), X = as.matrix(df.cov), 
                                      qs = seq(0.25, 0.75, by = 0.05), 
                                      q.fixed = 0.5, method = "exact")

# Merge
risk.overall_int_pfas <- risk.overall_int_pfas %>% mutate(outcome = "Internalizing") 
risk.overall_ext_pfas <- risk.overall_ext_pfas  %>% mutate(outcome = "Externalizing")
risk.overall_pfas <- rbind(risk.overall_int_pfas, risk.overall_ext_pfas)

## PBDE
df.exp <- mixture %>% dplyr::select(PBDE47.log, PBDE99.log, PBDE100.log)

# Internalizing
risk.overall_int_pbde <- OverallRiskSummaries(fit = fit_int_pbde, y = mixture$Internalizing, Z = as.matrix(df.exp), X = as.matrix(df.cov), 
                                      qs = seq(0.25, 0.75, by = 0.05), 
                                      q.fixed = 0.5, method = "exact")

# Externalizing

risk.overall_ext_pbde <- OverallRiskSummaries(fit = fit_ext_pbde, y = mixture$Externalizing, Z = as.matrix(df.exp), X = as.matrix(df.cov), 
                                      qs = seq(0.25, 0.75, by = 0.05), 
                                      q.fixed = 0.5, method = "exact")

# Merge
risk.overall_int_pbde <- risk.overall_int_pbde %>% mutate(outcome = "Internalizing") 
risk.overall_ext_pbde <- risk.overall_ext_pbde  %>% mutate(outcome = "Externalizing")
risk.overall_pbde <- rbind(risk.overall_int_pbde, risk.overall_ext_pbde)

# Plot

risk.overall_int_bkmr <- ggplot(risk.overall_int, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) + 
  geom_pointrange(size = 1.5, linewidth = 1.5, color = "grey25") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
  scale_y_continuous(limits = c(-0.4, 0.4), breaks = c(-0.4, -0.2, 0, 0.2, 0.4)) + 
  ggtitle("Total Sample
Internalizing") +
  theme_classic(base_family = "AppleGothic") +
  labs(x = "Quantile", y = "Beta (95 % Credible Interval)") +
  theme(text = element_text(size = 25, color = "grey25"),
        axis.text = element_text(color = "grey25", size = 25), 
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5))
risk.overall_ext_bkmr <- ggplot(risk.overall_ext, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) + 
  geom_pointrange(size = 1.5, linewidth = 1.5, color =  "grey55") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
  scale_x_continuous(limits = c(0.25, 0.75), breaks = c(0.30, 0.40, 0.50, 0.60, 0.70)) + 
  scale_y_continuous(limits = c(-0.4, 0.4), breaks = c(-0.4, -0.2, 0, 0.2, 0.4)) + 
  theme_classic() +
  labs(x = "Quantile", y = "Beta (95 % Credible Interval)", title = "Externalizing") +
  theme(text = element_text(size = 25, color = "grey25"),
        axis.text = element_text(color = "grey25", size = 25), 
        legend.title = element_blank(),
         plot.title = element_text(hjust = 0.5))
risk.overall_int_pfas_bkmr <- ggplot(risk.overall_int_pfas, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) + 
  geom_pointrange(size = 1.5, linewidth = 1.5, color = "grey25") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
  scale_y_continuous(limits = c(-0.4, 0.4), breaks = c(-0.4, -0.2, 0, 0.2, 0.4)) + 
  ggtitle("PFAS
     Internalizing") +
  theme_classic(base_family = "AppleGothic") +
  labs(x = "Quantile", y = "Beta (95 % Credible Interval)") +
  theme(text = element_text(size = 25, color = "grey25"),
        axis.text = element_text(color = "grey25", size = 25), 
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5))
risk.overall_ext_pfas_bkmr <- ggplot(risk.overall_ext_pfas, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) + 
  geom_pointrange(size = 1.5, linewidth = 1.5, color = "grey55") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
  scale_y_continuous(limits = c(-0.4, 0.4), breaks = c(-0.4, -0.2, 0, 0.2, 0.4)) + 
  theme_classic(base_family = "AppleGothic") +
  labs(x = "Quantile", y = "Beta (95 % Credible Interval)", title = "Externalizing") +
  theme(text = element_text(size = 25, color = "grey25"),
        axis.text = element_text(color = "grey25", size = 25), 
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5))
risk.overall_int_pbde_bkmr <- ggplot(risk.overall_int_pbde, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) + 
  geom_pointrange(size = 1.5, linewidth = 1.5, color = "grey25") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
  scale_y_continuous(limits = c(-0.4, 0.4), breaks = c(-0.4, -0.2, 0, 0.2, 0.4)) + 
  ggtitle("PBDE
    Internalizing") +
  theme_classic(base_family = "AppleGothic") +
  labs(x = "Quantile", y = "Beta (95 % Credible Interval)") +
  theme(text = element_text(size = 25, color = "grey25"),
        axis.text = element_text(color = "grey25", size = 25), 
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5))
risk.overall_ext_pbde_bkmr <- ggplot(risk.overall_ext_pbde, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) + 
  geom_pointrange(size = 1.5, linewidth = 1.5, color = "grey55") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
  scale_y_continuous(limits = c(-0.4, 0.4), breaks = c(-0.4, -0.2, 0, 0.2, 0.4)) + 
  theme_classic(base_family = "AppleGothic") +
  labs(x = "Quantile", y = "Beta (95 % Credible Interval)", title = "Externalizing") +
  theme(text = element_text(size = 25, color = "grey25"),
        axis.text = element_text(color = "grey25", size = 25), 
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5))
figure_4 <- plot_grid(risk.overall_int_bkmr, risk.overall_int_pfas_bkmr, risk.overall_int_pbde_bkmr, risk.overall_ext_bkmr,risk.overall_ext_pfas_bkmr, risk.overall_ext_pbde_bkmr, ncol = 3, row = 2)
figure_4

Figure S2

Univariate exposure-response functions indicating change in internalizing and externalizing behaviors resulting from individual natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in first trimester serum samples, while fixing remaining exposures in the mixture at the 50th percentile, estimated using Bayesian Kernal Machine Regression.

Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.

# Build labels
chemicals_label <- c("PFHXS", "PFOS", "PFOA", "PFNA", "BDE-47", "BDE-99", "BDE-100")
names(chemicals_label) <- c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")

# Internalizing
pred.resp.univar <- PredictorResponseUnivar(fit = fit_int)

int_univar <- ggplot(pred.resp.univar, aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) + 
  geom_smooth(stat = "identity", color  ="grey25", size = 2) + 
  facet_wrap(~ variable, label = labeller(variable = chemicals_label)) +
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  ylab("Internalizing") +
  xlab("Chemical") +
  theme_bw(base_family = "AppleGothic") +
  theme(text = element_text(size = 30, color = "grey25"),
      axis.text = element_text(color = "grey25", size = 30),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

# Externalizing
pred.resp.univar <- PredictorResponseUnivar(fit = fit_ext)


ext_univar <- ggplot(pred.resp.univar, aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) + 
  geom_smooth(stat = "identity", color = "grey25", size = 2) + 
  facet_wrap(~ variable, label = labeller(variable = chemicals_label)) +
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  xlab("Chemical") +
  ylab("Externalizing") +
  theme_bw(base_family = "AppleGothic") +
  theme(text = element_text(size = 30, color = "grey25"),
      axis.text = element_text(color = "grey25", size = 30),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

figure_s2 <- plot_grid(int_univar, ext_univar, ncol = 2, align = "h")
figure_s2

Figure S3

Bivariate exposure-response functions indicating change in internalizing and externalizing behaviors resulting from individual natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in first trimester serum samples, while fixing remaining exposures in the mixture at specific percentiles, estimated using Bayesian Kernal Machine Regression.

Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.

# Define mixture
df.exp <- mixture %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log)

# Internalizing

pred.resp.bivar <- PredictorResponseBivar(fit = fit_int, min.plot.dist = 1)
pred.resp.bivar.levels <- PredictorResponseBivarLevels(
  pred.resp.df = pred.resp.bivar, 
                                                       
  Z = as.matrix(df.exp), qs = c(0.25, 0.5, 0.75))

# Plot

# Create labels
chemicals_label <- c("PFHXS", "PFOS", "PFOA", "PFNA", "BDE-47", "BDE-99", "BDE-100")
names(chemicals_label) <- c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")

pred.resp.bivar.levels$variable1 <- factor(pred.resp.bivar.levels$variable1, levels = c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))
pred.resp.bivar.levels$variable2 <- factor(pred.resp.bivar.levels$variable2, levels = c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))
 

figure_s21 <- ggplot(pred.resp.bivar.levels, aes(z1, est)) + 
    geom_smooth(aes(col = quantile), stat = "identity") + 
    facet_grid(variable2 ~ variable1, label = labeller(variable1 = chemicals_label, variable2 = chemicals_label)) +
 scale_y_continuous(breaks = c(-0.4, 0, 0.4)) + 
  xlab("Chemical") +
  ylab("Internalizing") +
  scale_color_manual(values = c("steelblue3", "salmon3", "springgreen4")) +
  theme_bw(base_family = "AppleGothic") +
  theme(text = element_text(size = 30, color = "grey25"),
      axis.text = element_text(color = "grey25", size = 30),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
       legend.position = "bottom")
# Externalizing
pred.resp.bivar <- PredictorResponseBivar(fit = fit_ext, min.plot.dist = 1)
pred.resp.bivar.levels <- PredictorResponseBivarLevels(
  pred.resp.df = pred.resp.bivar, 
                                                       
  Z = as.matrix(df.exp), qs = c(0.25, 0.5, 0.75))

pred.resp.bivar.levels$variable1 <- factor(pred.resp.bivar.levels$variable1, levels = c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))
pred.resp.bivar.levels$variable2 <- factor(pred.resp.bivar.levels$variable2, levels = c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))

# Plot

figure_s22 <- ggplot(pred.resp.bivar.levels, aes(z1, est)) + 
    geom_smooth(aes(col = quantile), stat = "identity") + 
    facet_grid(variable2 ~ variable1, label = labeller(variable1 = chemicals_label, variable2 = chemicals_label)) +
  scale_y_continuous(breaks = c(-0.4, 0, 0.4)) + 
  xlab("Chemical") +
  ylab("Externalizing") +
  scale_color_manual(values = c("steelblue3", "salmon3", "springgreen4")) +
  theme_bw(base_family = "AppleGothic") +
  theme(text = element_text(size = 30, color = "grey25"),
      axis.text = element_text(color = "grey25", size = 30),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
       legend.position = "bottom")


figure_s3 <- plot_grid(figure_s21,figure_s22, col = 2, align = "h")
figure_s3

Self-Organizing Maps

Figure 4

Self-organizing maps (SOM) cluster star plot reflecting patterns of exposure to per- and polyfluoroalkyl substances (PFAS) and polybrominated diphenyl ethers (PBDEs), measured in serum samples obtained during the first trimester, and association between SOM clusters and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.

# Define mixture
df.exp <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log)

df.exp <- df.exp %>% rename(PFHXS = PFHXS.log, PFOS = PFOS.log, PFOA = PFOA.log, PFNA = PFNA.log, "BDE-47" = PBDE47.log, "BDE-99" = PBDE99.log, "BDE-100" = PBDE100.log)

# Explore clustering stats
# grp.str<-grpeval(df.exp, kmx=10)

# Initialize neurons/size of grid
grid <- somgrid(xdim = 3, ydim = 1, topo = "hexagonal")
# Train SOM model
set.seed(123)
som.mix <- som(X = as.matrix(df.exp), grid = grid, rlen = 500)

# Check if number of iterations is ok (should be stable) 
# plot(som.mix, type = "changes")

## Visualize

# N
# plot(som.mix, type = "counts")
# plot(som.mix, type = "mapping", keepMargins = TRUE)

# Codes/weights
par(cex = 8)
plot(som.mix, type = "codes", main = "SOM Clusters")

# Get unit classification
mix$unit_clusters <- som.mix$unit.classif

# get events (weight vector)
# getCodes(som.mix)

# Regression based on unit classifications

# Assign reference group 
mix$unit_clusters <- relevel(as.factor(mix$unit_clusters), ref = "1")

# Internalizing
som_int_reg <- lm(Internalizing ~ unit_clusters + `Mean Child Age`  + `Infant Sex` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use`, data = mix)

som_int_reg_df <- tidy(som_int_reg , conf.int = TRUE) %>% filter(term %in% c("unit_clusters2", "unit_clusters3")) %>% dplyr::select(term, estimate, conf.low, conf.high) %>% mutate(outcome = "Internalizing")


# Externalizing 
som_ext_reg <- lm(Externalizing ~ unit_clusters + `Mean Child Age`  + `Infant Sex` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use`, data = mix)

som_ext_reg_df <- tidy(som_ext_reg , conf.int = TRUE) %>% filter(term %in% c("unit_clusters2", "unit_clusters3")) %>% dplyr::select(term, estimate, conf.low, conf.high) %>% mutate (outcome = "Externalizing")


som_reg_df <- rbind(som_int_reg_df, som_ext_reg_df )
som_reg_df$outcome <- factor(som_reg_df$outcome, levels = c("Internalizing",  "Externalizing"))

figure_5 <- ggplot(som_reg_df, aes(x = estimate, y = term, color = outcome)) +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high), size = 4, linewidth = 4) +
  geom_vline(xintercept = 0, color = "grey25", linetype = "dashed", alpha = 0.5) +  
  facet_grid(~outcome, scales = "free") +
  xlab("Beta (95% Confidence Interval)") +
  ylab("Chemical") +
  scale_x_continuous(n.breaks=4) +
  scale_y_discrete(labels = c("Cluster 2", "Cluster 3")) +
  scale_color_manual(values = c("grey25", "grey55")) +
  coord_flip() +
  theme_classic(base_family = "AppleGothic") +
  theme(title = element_text(size = 50, color = "grey25"),
        axis.title.x = element_text(size = 50, color = "grey25"),
        axis.title.y = element_text(size = 50, color = "grey25"),
        axis.text.x = element_text(size = 50, color = "grey25", vjust = 0.6),
        axis.text.y = element_text(size = 50, color = "grey25"),
        strip.background = element_rect(fill="white"),
        strip.text = element_text(colour = 'grey25', face = "bold", size = 50), 
        legend.position = "none", 
        strip.background.x = element_blank(),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank())

figure_5

Table S7

Linear regression associations between self-organizing map (SOM) clusters and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.

Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.

# Table
table_s8 <- som_reg_df %>%
  mutate_if(is.numeric, round, digits = 2) %>%
  mutate(Beta_CI = paste0(estimate, " (", conf.low, ", ", conf.high, ")"))  %>%
 mutate(term = recode(term, 
                     "unit_clusters2" = "Cluster 2",
                     "unit_clusters3" = "Cluster 3")) %>%
  dplyr::select(term, outcome, Beta_CI) %>%
  pivot_wider(names_from = outcome, values_from = (Beta_CI)) %>%
  flextable() %>%
  add_body_row(values = c("" ,"Beta (95% CI)", "Beta (95% CI)")) %>%
  font(fontname = "Tahoma") %>%
  set_header_labels(term = "", sex = "Sex") %>%
  autofit() %>%
  theme_alafoli()

table_s8

	Internalizing	Externalizing
	Beta (95% CI)	Beta (95% CI)
Cluster 2	0.41 (-0.01, 0.83)	0.5 (0.09, 0.91)
Cluster 3	0.44 (0.08, 0.81)	0.54 (0.18, 0.89)